// /////////////////////////////////////////////////////////
// This file contains compatibility functions that allow
// FermiQCD to read and write in the Jim Simone data format
// Written by Massimo Di Pierro (2001)
// Requires MDP 1.3.4
//
// Please note that this is not te  native FermiQCD format
// To use the native format please use
//   generic_field::save()
//   generic_field::load()
// /////////////////////////////////////////////////////////
#ifndef __JIM_IO_H__
#define __JIM_IO_H__

typedef std::map<string,string> MetaHash;

template<class T>
struct jim_complex {
  T re;
  T im;
};

// /////////////////////////////////////////////////////////
// Macro and structures used by the JIM data format
// /////////////////////////////////////////////////////////

#define JIM_IO_UNI_MAGIC 0x71626434
typedef long JIM_IO_type32;

typedef struct
{
  JIM_IO_type32 magic_number;
  JIM_IO_type32 gmtime;
  JIM_IO_type32 size_of_element;
  JIM_IO_type32 elements_per_site;
  JIM_IO_type32 dim[4];         /*       { X=0, Y, Z, T}      */
  JIM_IO_type32 site_order;     /* 0 for [t][z][y][x] order   */
} JIM_header;


void JIM_read_metadata ( const string& basefile, MetaHash& meta,
                         const int IOprocessRank = 0,
                         const unsigned maxlen = 2048 )
{
  if ( ME == IOprocessRank )
    {
      string file = basefile;
      file += ".info";

      ifstream mstream;
      mstream.open ( file.c_str ( ) );

      const string white ( " \t\n" );
      while ( 1 )
        {
          char buf [maxlen];

          mstream.getline ( buf, maxlen );

          if ( mstream.eof ( ) )
            {
              break;
            }

          string line ( buf );

          int indx = line.find_first_of ( white );

          string key = line.substr ( 0, indx );

          indx = line.find_first_not_of ( white, indx );

          string value = line.substr ( indx, line.length ( ) );

          meta[key] = value;
        }

      mstream.close ( );
    }
}

void JIM_print_metadata ( ostream& os, const MetaHash& meta )
{
  MetaHash::const_iterator j;
  for ( j = meta.begin ( ); j != meta.end ( ); ++j )
    {
      os << j -> first << "\t\t" << j -> second << "\n";
    }
}

template<class T>
string stringify ( const T value, const char* format, const int maxlen = 256 )
{
  char buf [maxlen];

  sprintf ( buf, format, value );

  return string ( buf );
}

void get_header_meta ( const JIM_header& header, MetaHash& meta )
{
  meta ["header.magic-number"] =
    stringify ( header.magic_number, "0x%08x" );
 
  meta ["header.timestamp"] =
    stringify ( ctime ( static_cast<const time_t *> ( &header.gmtime ) ),
                "%s" );

  meta ["header.size-of-element"] =
    stringify ( header.size_of_element, "%d" );

  meta ["header.elements-per-site"] =
    stringify ( header.elements_per_site, "%d" );

  meta ["header.dim[0]"] =
    stringify ( header.dim[0], "%d" );

  meta ["header.dim[1]"] =
    stringify ( header.dim[1], "%d" );

  meta ["header.dim[2]"] =
    stringify ( header.dim[2], "%d" );

  meta ["header.dim[3]"] =
    stringify ( header.dim[3], "%d" );

  meta ["header.site-order"] =
    stringify ( header.site_order, "%d" );
}

void JIM_write_metadata ( const string filename,
                          const JIM_header& header, const MetaHash& user,
                          const int IOprocessRank = 0 )
{
  // Write Metadata
  if ( ME == IOprocessRank )
    {
      MetaHash Mlist;

      get_header_meta ( header, Mlist );

      string metafile = filename;
      metafile += ".info";
      ofstream smeta;
      smeta.open ( metafile.c_str ( ) );
      if ( smeta.fail ( ) )
        {
	  std::cout << "write_metadata: open failed" << "\n";
        }

      MetaHash::const_iterator j;
      for ( j = Mlist.begin ( ); j != Mlist.end ( ); ++j )
        {
          smeta << j -> first << "\t\t" << j -> second << "\n";
        }

      // Write user metadata
      for ( j = user.begin ( ); j != user.end ( ); ++j )
        {
          smeta << j -> first << "\t\t" << j -> second << "\n";
        }

      smeta.close ( );
    }
}

void JIM_append_metadata ( const string filename,
                           const MetaHash& user,
                           const int IOprocessRank = 0 )
{
  if ( ME == IOprocessRank )
    {
      string metafile = filename;
      metafile += ".info";
      ofstream smeta;

      smeta.open ( metafile.c_str ( ), ios::out | ios::app );
      if ( smeta.fail ( ) )
        {
	  std::cout << "write_metadata: open+append failed" << "\n";
        }

      MetaHash::const_iterator j;
      // Write user metadata
      for ( j = user.begin ( ); j != user.end ( ); ++j )
        {
          smeta << j -> first << "\t\t" << j -> second << "\n";
        }

      smeta.close ( );
    }
}

void JIM_print_header ( const JIM_header& header,
                        const int IOprocessRank = 0 )
{
  if ( ME == IOprocessRank )
    {

      printf ( "header.magic-number\t0x%8lx\n",   header.magic_number );
      printf ( "header.gmtime\t%s",
               ctime ( &header.gmtime ) );
      printf ( "header.size-of-element\t%ld\n",   header.size_of_element );
      printf ( "header.elements_per-site\t%ld\n", header.elements_per_site );
      printf ( "header.dim[0]\t%ld\n",            header.dim[0] );
      printf ( "header.dim[1]\t%ld\n",            header.dim[1] );
      printf ( "header.dim[2]\t%ld\n",            header.dim[2] );
      printf ( "header.dim[3]\t%ld\n",            header.dim[3] );
      printf ( "header.site-order\t%ld\n",        header.site_order );
    }
}

bool JIM_header_mismatch ( const JIM_header& header0,
                           const JIM_header& header1 )
{
  return
    header0.magic_number      != header1.magic_number ||
    header0.size_of_element   != header1.size_of_element ||
    header0.elements_per_site != header1.elements_per_site ||
    header0.dim[0]            != header1.dim[0] ||
    header0.dim[1]            != header1.dim[1] ||
    header0.dim[2]            != header1.dim[2] ||
    header0.dim[3]            != header1.dim[3] ||
    header0.site_order        != header1.site_order;
}

void JIM_reverse_header_endianess ( JIM_header& header )
{
  switch_endianess_byte4 ( header.magic_number );

  if ( header.magic_number != JIM_IO_UNI_MAGIC )
    {
      std::cout << "wrong header magic after byte swap" << "\n";
      exit ( 1 );
    }

  switch_endianess_byte4 ( header.size_of_element );
  switch_endianess_byte4 ( header.elements_per_site );
  switch_endianess_byte4 ( header.dim [0] );
  switch_endianess_byte4 ( header.dim [1] );
  switch_endianess_byte4 ( header.dim [2] );
  switch_endianess_byte4 ( header.dim [3] );
  switch_endianess_byte4 ( header.site_order );
}

// /////////////////////////////////////////////////////////
// This function is used by mdp_field::save()/load()
// to sort the sites before/after I/O to conform to
// JIM data format
// /////////////////////////////////////////////////////////

// internal q [t][x][y][z]
// external q [t][z][y][x]
long JIM_sort ( mdp_lattice &lattice, long idx )
{
  int x [4], tmp;

  lattice.global_coordinate ( idx, x );

  /* swap X and Z */
  tmp         = x [SPACE_X];
  x [SPACE_X] = x [SPACE_Z];
  x [SPACE_Z] = tmp;
  return lattice.global_coordinate ( x );
}

// /////////////////////////////////////////////////////////
// Parallel I/O compatible with JIM file format
// Field: gauge_field
// /////////////////////////////////////////////////////////

template<class T>
void JIM_save ( gauge_field &U, char *filename, MetaHash& user )
{
  JIM_header header;
  int i,mu,nu,a,b, nc=U.nc;
  mdp_field<jim_complex<T> > new_u(U.lattice(),4*nc*nc);
  int IOprocessRank = 0;
  long BufSize = 512;

  site x(U.lattice());

  if ( U.ndim != 4 ) error("JIM IO only for ndim=4");

  header.magic_number=JIM_IO_UNI_MAGIC;
  header.size_of_element=  sizeof(T);
  header.elements_per_site=8*nc*nc;
  header.dim[0]=U.lattice().size(SPACE_X);
  header.dim[1]=U.lattice().size(SPACE_Y);
  header.dim[2]=U.lattice().size(SPACE_Z);
  header.dim[3]=U.lattice().size(TIME);
  time_t time_and_date;
  header.gmtime = time ( &time_and_date );
  header.site_order = 0;

  // direction: T X Y Z
  //  internal: 0 1 2 3
  //  file    : 3 0 1 2
  int mapdir [4] = { 3, 0, 1, 2 };

  forallsites ( x )
    {
      for ( mu = 0; mu < 4; mu++ )
        {
          nu = mapdir [mu];

          for ( a = 0; a < nc; a++ )
            for ( b = 0; b < nc; b++ )
              {
                i = ( ( nu * nc ) + a ) * nc + b;
                new_u ( x, i ).re = U ( x, mu, a, b ).re;
                new_u ( x, i ).im = U ( x, mu, a, b ).im;
	      }
        }   
    }

  if ( ME == IOprocessRank && shutup != FALSE )
    printf("Saving a gauge field in JIM  format.\n");

  new_u.save ( filename, IOprocessRank, BufSize,
	       (char*) &header,
               sizeof(JIM_header), JIM_sort );

  JIM_write_metadata ( string (filename), header, user, IOprocessRank );
}

// Used for wavefunction
template<class T>
void JIM_load ( mdp_field<mdp_complex>& wf, const char* filename,
                MetaHash& meta )
{
  JIM_header header;
  int IOprocessRank = 0;
  long BufSize = 512;
  mdp_field<jim_complex<T> > _wf ( wf.lattice ( ) ); // precision T
  site x ( _wf.lattice ( ) );

  if ( _wf.lattice ( ).ndim != 4 )
    error ( "JIM IO only for ndim == 4" );

  if ( ME == IOprocessRank && shutup != FALSE )
    printf ( "Loading a mdp_field<mdp_complex> field in JIM  format.\n" );

  // load() does not switch endianess automatically if not NATIVE mode
  _wf.load ( (char *) filename, IOprocessRank, BufSize,
             (char *) &header,
             sizeof (JIM_header), JIM_sort );

  // valid header only on node IOprocessRank
  if ( ME == IOprocessRank )
    {
      if ( header.magic_number != JIM_IO_UNI_MAGIC )
        {
	  std::cout << "byte swapping file" << "\n";
          JIM_reverse_header_endianess ( header );
          _wf.switch_endianess_4bytes ( );
        }

      //DEBUG
      //JIM_print_header ( header, IOprocessRank );

      if ( header.size_of_element != sizeof (T) )
        {
          error ( "IO: Wrong header.size-of-element!" );
        }
      if ( header.elements_per_site != 2 )
        {
          error ( "IO: Wrong header.elements-per-site!" );
        }
      if ( header.dim [0] != _wf.lattice ( ).size (SPACE_X) )
        {
          error ( "IO: Wrong header.dim[0]" );
        }
      if ( header.dim [1] != _wf.lattice ( ).size (SPACE_Y) )
        {
          error("IO: Wrong header.dim[1]");
        }
      if ( header.dim [2] != _wf.lattice ( ).size (SPACE_Z) )
        {
          error ( "IO: Wrong header.dim[2]" );
        }
      if ( header.dim [3] != _wf.lattice ( ).size (TIME) )
        {
          error ( "IO: Wrong header.dim[3]" );
        }
      if ( header.site_order != 0 )
        {
          error ( "IO: Unknown header.site-order" );
        }
    }

  // copy to convert from T to default precision
  forallsitesandcopies ( x )
    {
      wf ( x ).re = _wf ( x ).re;
      wf ( x ).im = _wf ( x ).im;
    }

  // get metadata
  JIM_read_metadata ( string (filename), meta, IOprocessRank );
}

template<class T>
void JIM_load ( gauge_field &U, char *filename, MetaHash& gaugeMeta )
{
  JIM_header header;
  int IOprocessRank = 0;
  long BufSize = 512;
  mdp_field<jim_complex<T> > new_u(U.lattice(),4*U.nc*U.nc);
  site x(U.lattice());
  int i,mu,nu,a,b, nc=U.nc;

  if(U.ndim!=4) error("JIM IO only for ndim=4");
  if ( ME == 0 && shutup != FALSE )
    printf("Loading a gauge field in JIM  format.\n");

  // load() does not switch endianess automatically if not NATIVE mode
  new_u.load ( filename, IOprocessRank, BufSize,
               (char *) &header,
               sizeof(JIM_header), JIM_sort);

  // valid header only on node IOprocessRank
  if ( ME == IOprocessRank )
    {
      if ( header.magic_number != JIM_IO_UNI_MAGIC )
        {
          JIM_reverse_header_endianess ( header );
          new_u.switch_endianess_4bytes ( );
        }

      if(header.size_of_element!=  sizeof(T))
        error("IO: Wrong size_of_element!");
      if(header.elements_per_site!=8*nc*nc)
        error("IO: Wrong elements_per_site!");
      if(header.dim[0]!=U.lattice().size(SPACE_X))
        error("IO: Wrong ndim[0]");
      if(header.dim[1]!=U.lattice().size(SPACE_Y))
        error("IO: Wrong ndim[1]");
      if(header.dim[2]!=U.lattice().size(SPACE_Z))
        error("IO: Wrong ndim[2]");
      if(header.dim[3]!=U.lattice().size(TIME))
        error("IO: Wrong ndim[3]");
      if(header.site_order!=0)
        error("IO: Unknown site order");
    }

  // direction: T X Y Z
  //  internal: 0 1 2 3
  //  file    : 3 0 1 2
  int mapdir [4] = { 3, 0, 1, 2 };

  forallsitesandcopies(x) {
    for(mu=0; mu<4; mu++) {

      nu = mapdir [mu];

      for ( a = 0; a < nc; a++ )
        for ( b = 0; b < nc; b++ )
          {
            i = ( ( nu * nc ) + a ) * nc + b;
            U ( x, mu, a, b ).re = new_u ( x, i ).re;
            U ( x, mu, a, b ).im = new_u ( x, i ).im;
	  }
    }
  }

  // get metadata
  JIM_read_metadata ( string (filename), gaugeMeta, IOprocessRank );
}

// /////////////////////////////////////////////////////////
// Parallel I/O compatible with Jim Simone file format
// Field: fermi_field
// /////////////////////////////////////////////////////////

// TODO: User meta, checksums, record list
template<class T>
void JIM_save ( fermi_field &psi, char *filename, MetaHash& user,
                int spin = 0, int color = 0 )
{
  int IOprocessRank = 0;
  long BufSize = 512;
  JIM_header header;
  int i,a,b, nc=psi.nc;
  mdp_field<jim_complex<T> > new_psi(psi.lattice(),4*nc);
  site x(psi.lattice());

  if(psi.lattice().ndim!=4) error("JIM IO only for ndim=4");

  header.magic_number=JIM_IO_UNI_MAGIC;
  header.size_of_element=  sizeof(T);
  header.elements_per_site=8*nc;
  header.dim[0]=psi.lattice().size(SPACE_X);
  header.dim[1]=psi.lattice().size(SPACE_Y);
  header.dim[2]=psi.lattice().size(SPACE_Z);
  header.dim[3]=psi.lattice().size(TIME);
  time_t time_and_date;
  header.gmtime=time(&time_and_date);
  header.site_order=0;

  forallsites(x) {
    for(a=0; a<4; a++) {
      for(b=0; b<nc; b++) {
        // JIM mapping!!
        i=a*nc+b;
        new_psi ( x, i ).re = psi ( x, a, b ).re;
        new_psi ( x, i ).im = psi ( x, a, b ).im;
      }
    }
  }
  if ( ME == 0 && shutup != FALSE )
    printf("Saving a fermi field in JIM  format.\n");

  if ( spin == 0 && color == 0 )
    {
      new_psi.save ( filename, IOprocessRank, BufSize,
                     ( char* ) &header,
                     sizeof(JIM_header), JIM_sort );

      // Write Metadata
      JIM_write_metadata ( string(filename), header, user, IOprocessRank );
    }
  else
    {
      if ( ME == IOprocessRank )
        {
          JIM_header checkheader;
          bool writehdr = false;
          FILE *fp = fopen ( filename, "rb+" );

          // check for valid header
          if ( 1 != fread ( &checkheader, 1, sizeof(JIM_header), fp ) )
            {
              writehdr = true;
            }
          else
            {
              if ( JIM_header_mismatch ( checkheader, header ) )
                {
		  std::cout << "header mismatch upon append!" << "\n";
                  exit ( 1 );
                }
            }

          if ( writehdr )
            {
              fseek ( fp, 0, SEEK_SET );
              fwrite ( &header, 1, sizeof(JIM_header), fp );
            }
          fclose(fp);
        }

      // write data
      new_psi.save ( filename, IOprocessRank, BufSize, "NOHEADER",
                     sizeof(JIM_header) + new_psi.field_size ( )
                     * ( spin * nc + color ),
                     JIM_sort, "r+b");
    }

  if ( ME == IOprocessRank )
    {
      // append meta data
      string metafile = filename;
      metafile += ".info";
      FILE *meta_fp = fopen ( metafile.c_str ( ), "a" );

      // record number
      fprintf ( meta_fp, "record[%d].spin-color\t\t\t%d,%d\n",
                3 * spin + color, spin, color );

      fclose ( meta_fp );
    }
}


template<class T>
void JIM_load ( fermi_field &psi, char *filename,
		MetaHash& QuarkMeta, int spin = 0, int color = 0 )
{
  JIM_header header;
  int i,a,b, nc=psi.nc;
  mdp_field<jim_complex<T> > new_psi(psi.lattice(),4*nc);
  site x(psi.lattice());
  int IOprocessRank = 0;
  long BufSize = 512;

  if(psi.lattice().ndim!=4)
    error("JIM IO only for ndim=4");

  if(ME==0 && shutup != FALSE)
    printf("Loading a fermi field in JIM  format.\n");

  if ( ME == IOprocessRank )
    {
      // valid header only on node IOprocessRank

      FILE *fp = fopen ( filename, "rb" );
      fread ( &header, 1, sizeof(JIM_header), fp );
      fclose ( fp );

      if ( header.magic_number != JIM_IO_UNI_MAGIC )
        {
          JIM_reverse_header_endianess ( header );
          new_psi.switch_endianess_4bytes ( );
        }

      if(header.size_of_element!=  sizeof(T))
        error("Wrong size!");
      if(header.elements_per_site!=8*nc)
	error("Wrong nc!");
      if(header.dim[0]!=psi.lattice().size(SPACE_X))
	error("Wrong ndim[0]");
      if(header.dim[1]!=psi.lattice().size(SPACE_Y))
	error("Wrong ndim[1]");
      if(header.dim[2]!=psi.lattice().size(SPACE_Z))
	error("Wrong ndim[2]");
      if(header.dim[3]!=psi.lattice().size(TIME))
	error("Wrong ndim[3]");
         
      if(header.site_order!=0)
	error("Unkown site order");

    }

  // load() does not switch endianess automatically if not NATIVE mode
  if ( spin == 0 && color == 0 )
    {
      new_psi.load ( filename, IOprocessRank, BufSize,
		     (char*) &header,
		     sizeof(JIM_header), JIM_sort );
    }
  else
    {
      new_psi.load ( filename, IOprocessRank, BufSize,
                     "NOHEADER",
                     sizeof(JIM_header) +
                     new_psi.field_size ( ) * ( spin * nc + color ),
                     JIM_sort );
    }

  forallsitesandcopies(x) {
    for(a=0; a<4; a++)
      for(b=0; b<nc; b++) {
        // JIM mapping!!
        i=a*nc+b;
        psi ( x, a, b ).re = new_psi ( x, i ).re;
        psi ( x, a, b ).im = new_psi ( x, i ).im;
      }
  }

  // get metadata
  JIM_read_metadata ( string (filename), QuarkMeta, IOprocessRank );
}

// /////////////////////////////////////////////////////////
// Parallel I/O compatible with Jim Simone file format
// Field: fermi_propagator
// /////////////////////////////////////////////////////////


template<class T>
void JIM_save ( fermi_propagator &S, char *filename, MetaHash& user )
{
  JIM_header header;
  int i,a,b,mu,nu,nc=S.nc;
  mdp_field<jim_complex<T> > new_S(S.lattice(),16*nc*nc);
  site x(S.lattice());
  int IOprocessRank = 0;
  long BufSize = 512;

  if(S.lattice().ndim!=4) error("JIM IO only for ndim=4");

  header.magic_number=JIM_IO_UNI_MAGIC;
  header.size_of_element=  sizeof(T);
  header.elements_per_site=32*nc*nc;
  header.dim[0]=S.lattice().size(SPACE_X);
  header.dim[1]=S.lattice().size(SPACE_Y);
  header.dim[2]=S.lattice().size(SPACE_Z);
  header.dim[3]=S.lattice().size(TIME);
  time_t time_and_date;
  header.gmtime=time(&time_and_date);
  header.site_order=0;

  forallsites(x) {
    for(mu=0; mu<4; mu++)
      for(nu=0; nu<4; nu++)
        for(a=0; a<nc; a++)
          for(b=0; b<nc; b++) {
            // JIM mapping!!
            i=(((mu*nc)+a)*4+nu)*nc+b;
            new_S ( x, i ).re = S ( x, mu, nu, a, b ).re;
            new_S ( x, i ).im = S ( x, mu, nu, a, b ).im;
          }
  }
  if(ME==0 && shutup != FALSE)
    printf("Saving a fermi propagator in JIM  format.\n");

  new_S.save ( filename, IOprocessRank, BufSize,
	       (char*) &header, sizeof(JIM_header), JIM_sort );

  // Write Metadata
  if ( ME == IOprocessRank )
    {
      MetaHash Mlist;

      get_header_meta ( header, Mlist );

      string metafile = filename;
      metafile += ".info";
      FILE *meta_fp = fopen ( metafile.c_str ( ), "w" );
      MetaHash::iterator j;
      for ( j = Mlist.begin ( ); j != Mlist.end ( ); ++j )
        {
          fprintf ( meta_fp, "%s\t\t\t%s\n",
                    j -> first.c_str ( ),
                    j -> second.c_str ( ) );
        }

      // write matrix ordering
      fprintf ( meta_fp,
                "quark.matrix-order\t\t\tmdp_complex<%s> Q[src-spin][src-color][spin][color]\n",
                sizeof(mdp_complex) == 8 ? "float" : "double" );

      // Write user metadata
      for ( j = user.begin ( ); j != user.end ( ); ++j )
        {
          fprintf ( meta_fp, "%s\t\t\t%s\n",
                    j -> first.c_str ( ),
                    j -> second.c_str ( ) );
        }


      // TODO: checksum by timeslice ?

      fclose ( meta_fp );
    }
}


template<class T>
void JIM_load ( fermi_propagator &S, char *filename, MetaHash& meta )
{
  JIM_header header;
  int i,a,b,mu,nu,nc=S.nc;
  mdp_field<jim_complex<T> > new_S(S.lattice(),16*nc*nc);
  site x(S.lattice());
  int IOprocessRank = 0;
  long BufSize = 512;

  if(S.lattice().ndim!=4)
    error("JIM IO only for ndim=4");

  if(ME==0 && shutup != FALSE)
    printf("Loading a fermi propagator in JIM  format.\n");

  // load() does not switch endianess automatically if not NATIVE mode
  new_S.load ( filename, IOprocessRank, BufSize, (char*) &header,
               sizeof(JIM_header), JIM_sort );

  // valid header only on node IOprocessRank
  if ( ME == IOprocessRank )
    {
      if ( header.magic_number != JIM_IO_UNI_MAGIC )
        {
          JIM_reverse_header_endianess ( header );
          new_S.switch_endianess_4bytes ( );
        }

      if(header.size_of_element!=  sizeof(T))
        error("Wrong size!");
      if(header.elements_per_site!=32*nc*nc)
        error("Wrong nc!");
      if(header.dim[0]!=S.lattice().size(SPACE_X))
        error("Wrong ndim[0]");
      if(header.dim[1]!=S.lattice().size(SPACE_Y))
        error("Wrong ndim[1]");
      if(header.dim[2]!=S.lattice().size(SPACE_Z))
        error("Wrong ndim[2]");
      if(header.dim[3]!=S.lattice().size(TIME))
        error("Wrong ndim[3]");
      if(header.site_order!=0)
        error("Unkown site order");
    }

  forallsitesandcopies(x) {
    for(mu=0; mu<4; mu++)
      for(nu=0; nu<4; nu++)
        for(a=0; a<nc; a++)
          for(b=0; b<nc; b++) {
            // JIM mapping!!
            i=(((mu*nc)+a)*4+nu)*nc+b;
            S ( x, mu, nu, a, b ).re = new_S ( x, i ).re;
            S ( x, mu, nu, a, b ).im = new_S ( x, i ).im;
          }
  }

  // get metadata
  JIM_read_metadata ( string (filename), meta, IOprocessRank );
}

#ifdef RUN_TESTS
// //////////////////////////////////////////////////////
// Test program for  JIM_IO libs
// programs requires:
//    MDP_Lib2, MDP_MPI, MDP_Gauge, MDP_Fermi, JIM_IO
// //////////////////////////////////////////////////////

int main_test_JIM_IO(int argc, char **argv) {
  mpi.open_wormholes(argc,argv);

  int i=0,a,b,mu,nu, box[]={6,6,6,6};
  mdp_lattice      spacetime(4, box);
  gauge_field      U(spacetime, 3);
  fermi_field      psi(spacetime,3);
  fermi_field      chi(spacetime,3);
  fermi_propagator S(spacetime,3);
  site             x(spacetime);

  i=0;
  forallsites(x) {
    for(mu=0; mu<4; mu++)
      for(a=0; a<U.nc; a++)
        for(b=0; b<U.nc; b++) {
          U(x,mu,a,b)=mdp_complex(i,0);
          i++;
        }
  }
  i=0;
  forallsites(x) {
    for(mu=0; mu<4; mu++)
      for(a=0; a<U.nc; a++) {
        psi(x,mu,a)=mdp_complex(i,0);
        chi(x,mu,a)=mdp_complex(0,i);
        i++;
      }
  }
  i=0;
  forallsites(x) {
    for(mu=0; mu<4; mu++)
      for(nu=0; nu<4; nu++)
        for(a=0; a<U.nc; a++)
          for(b=0; b<U.nc; b++) {
            S(x,mu,nu,a,b)=mdp_complex(i,0);
            i++;
          }
  }

  JIM_save<float>(U, "test.dat");
  JIM_load<float>(U, "test.dat");
  JIM_save<float>(psi, "test.dat");
  JIM_save<float>(chi, "test.dat", 2, 1);
  JIM_load<float>(psi, "test.dat");
  JIM_load<float>(chi, "test.dat", 2, 1);
  JIM_save<float>(S, "test.dat");
  JIM_load<float>(S, "test.dat");

  i=0;
  forallsites(x) {
    for(mu=0; mu<4; mu++)
      for(a=0; a<U.nc; a++)
        for(b=0; b<U.nc; b++) {
          if(U(x,mu,a,b)!=mdp_complex(i,0)) error("U");
          i++;   
        }
  }
  printf("JIM I/O for gauge_field is ok\n");

  i=0;
  forallsites(x) {
    for(mu=0; mu<4; mu++)
      for(a=0; a<U.nc; a++) {
        if(psi(x,mu,a)!=mdp_complex(i,0)) error("psi");
        if(chi(x,mu,a)!=mdp_complex(0,i)) error("chi");
        i++;
      }
  }
  printf("JIM I/O for fermi_field is ok\n");

  i=0;
  forallsites(x) {
    for(mu=0; mu<4; mu++)
      for(nu=0; nu<4; nu++)
        for(a=0; a<U.nc; a++)
          for(b=0; b<U.nc; b++) {
            if(S(x,mu,nu,a,b)!=mdp_complex(i,0)) error("S");
            i++;
          }
  }
  printf("JIM I/O for fermi_propagator is ok\n");

  mpi.close_wormholes();
  return 0;
}
#endif

#endif

